/************************************************************************
 * MechSys - Open Library for Mechanical Systems                        *
 * Copyright (C) 2013 Sergio Torres                                     *
 *                                                                      *
 * This program is free software: you can redistribute it and/or modify *
 * it under the terms of the GNU General Public License as published by *
 * the Free Software Foundation, either version 3 of the License, or    *
 * any later version.                                                   *
 *                                                                      *
 * This program is distributed in the hope that it will be useful,      *
 * but WITHOUT ANY WARRANTY; without even the implied warranty of       *
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the         *
 * GNU General Public License for more details.                         *
 *                                                                      *
 * You should have received a copy of the GNU General Public License    *
 * along with this program. If not, see <http://www.gnu.org/licenses/>  *
 ************************************************************************/

// STL
#include <iostream>
#include <fstream>
#include <cmath>
#include <cstdlib> // for srand, rand
#include <ctime>   // for time

// MechSys
#include <mechsys/util/fatal.h>
#include <mechsys/util/array.h>
#include <mechsys/util/util.h>
#include <mechsys/linalg/matvec.h>
#include <mechsys/mesh/structured.h>
#include <mechsys/mesh/unstructured.h>
#include <mechsys/dem/domain.h>

using std::cout;
using std::endl;
using std::ifstream;
using std::ofstream;
using std::setw;
using std::min;
using std::max;

double determ(double cd[9])
{
	double det;
	det=cd[0]*cd[4]*cd[8]+cd[1]*cd[5]*cd[6]+cd[2]*cd[3]*cd[7];
	det=det-cd[2]*cd[4]*cd[6]-cd[1]*cd[3]*cd[8]-cd[0]*cd[5]*cd[7];
	return(det);
}


void cross(double cd[8],double cro[2])
{
	double bb[2];
	
	bb[0] = cd[0]*cd[5]-cd[4]*cd[1];
	bb[1] = cd[2]*cd[7]-cd[6]*cd[3];
	
	double det=(cd[0] - cd[1])*(cd[6] - cd[7])-(cd[4] - cd[5])*(cd[2] - cd[3]);
	
	cro[0] = bb[0]*(cd[2] - cd[3])-bb[1]*(cd[0] - cd[1]);
	cro[1] = bb[0]*(cd[6] - cd[7])-bb[1]*(cd[4] - cd[5]);
	cro[0] = cro[0]/det;
	cro[1] = cro[1]/det;
	

}

int main(int argc, char **argv) try
{
    if (argc<2) throw new Fatal("This program must be called with one argument: the name of the data input file without the '.inp' suffix.\nExample:\t %s filekey\n",argv[0]);
    String filekey  (argv[1]);
    String filename (filekey+".inp");
    ifstream infile(filename.CStr());



    /*Preprocessing the basic quantities*/
    /*Pi*/
    double Pi=3.14159265358979;
    /*dimensions - these need to remain as they are*/
    double Lx=5.0 , Ly=2.0 , Lz=0.2;
    /*dummy dimensions - these can be tampered with*/
    double Lx0=4. , Ly0=2.;
    /*Total number of elements - approximate*/
    size_t Ntot=10;
    /*scaling: dxdy=n means the block dimension x is n times the dimension y*/
    double dxdy=2.;
    /*mesh dimensions*/
    size_t Nx=2  , Ny=4  , Nz=1;
    /*element dimensions*/
    double lx=5.0 , ly=2.;
    /*inclination angles*/
    double angx=0.0 , angy=0.0;
    /*spatial dispersion*/
    double lsx=5.0 , lsy=2.;
    
    /*reading from file*/
    infile >> Lx;          infile.ignore(200,'\n');
    infile >> Ly;          infile.ignore(200,'\n');
    infile >> Lz;          infile.ignore(200,'\n');
    infile >> Ntot;        infile.ignore(200,'\n');
    infile >> dxdy;        infile.ignore(200,'\n');
    infile >> Nz;          infile.ignore(200,'\n'); 
    infile >> angx;        infile.ignore(200,'\n');
    infile >> angy;        infile.ignore(200,'\n');  
    infile >> lsx;         infile.ignore(200,'\n');
    infile >> lsy;         infile.ignore(200,'\n');  

    cout << angx<<" , "<<angy<<"\n";
    if ((angx+angy)>89) throw new Fatal("This selection of angles is not acceptable");
    if ((angx+angy)<-89) throw new Fatal("This selection of angles is not acceptable");

    /*converting angles to radians*/
    angx = angx*Pi/double(180.);
    angy = angy*Pi/double(180.);
  
    Ny=(int)round(sqrt((double)Ntot*Ly/Lx*dxdy));
    
    ly=Ly/(double)Ny;
    lx=ly*dxdy;

    Nx=(int)(Lx/lx);


    cout << "Nx=" << Nx << ", Ny=" << Ny <<"\n";
    cout << "lx=" << lx << ", ly=" << ly <<"\n";
    cout << "Lx0=" << Lx0 << ", Ly0=" << Ly0 <<"\n";
    cout << "Lx=" << Lx << ", Ly=" << Ly <<"\n";


    Vec3_t       cell0;       ///< X, Y, and Z coordinates
    Vec3_t       cell1;       ///< X, Y, and Z coordinates
    Vec3_t       cell2;       ///< X, Y, and Z coordinates
    Vec3_t       cell3;       ///< X, Y, and Z coordinates
    Vec3_t       L1;          ///< X, Y, and Z coordinates
    Vec3_t       L2;          ///< X, Y, and Z coordinates
    Vec3_t       L3;          ///< X, Y, and Z coordinates
    Vec3_t       L4;          ///< X, Y, and Z coordinates
    Vec3_t       c;           ///< X, Y, and Z coordinates
    double       a1,a2,b1,b2;

    cell0 = 0.0, 0.0, 0.0;
    cell1 = -ly*sin(angy), ly*cos(angy), 0.0;
    cell2 = -ly*sin(angy)+lx*cos(angx), ly*cos(angy)-lx*sin(angx), 0.0;
    cell3 = lx*cos(angx), -lx*sin(angx), 0.0;


    /*Evalating the mesh dimensions*/

    /*Setting the total length*/
    /*first block vertex*/
    a1=0.0;
    a2=0.0;
    b1=0.0;
    b2=0.0;
    /*second block vertex*/
    a1= min(a1,(sin(angx)*Lx));
    a2= max(a2,(sin(angx)*Lx));
    b1= min(b1,(cos(angy)*Lx));
    b2= max(b2,(cos(angy)*Lx));
    /*third block vertex*/
    a1= min(a1,(sin(angx)*Lx+cos(angx)*Ly));
    a2= max(a2,(sin(angx)*Lx+cos(angx)*Ly));
    b1= min(b1,(cos(angy)*Lx+sin(angy)*Ly));
    b2= max(b2,(cos(angy)*Lx+sin(angy)*Ly));
    /*fourth block vertex*/
    a1= min(a1,(cos(angx)*Ly));
    a2= max(a2,(cos(angx)*Ly));
    b1= min(b1,(sin(angy)*Ly));
    b2= max(b2,(sin(angy)*Ly));

    Nx=(b2-b1)/(cos(angx+angy)*lx);
    Ny=(a2-a1)/(cos(angx+angy)*ly);
    c[0]=(b2-b1)/(cos(angx+angy)*lx)-(double)Nx;
    c[1]=(a2-a1)/(cos(angx+angy)*ly)-(double)Ny;
    Nx++;
    Ny++;

    c[0]=(1.0-c[0])*(cos(angx+angy)*lx)/2.0;
    c[1]=(1.0-c[1])*(cos(angx+angy)*ly)/2.0;

    a1 = a1 - c[1];
    a2 = a2 + c[1];
    b1 = b1 - c[0];
    b2 = b2 + c[0];

    /*finding the intersection points*/
    L1[0]=-( sin(angy)*a1 -cos(angx)*b1)/cos(angx+angy);
    L1[1]=-(-cos(angy)*a1 +sin(angx)*b1)/cos(angx+angy);

    L2[0]=-( sin(angy)*a2 -cos(angx)*b1)/cos(angx+angy);
    L2[1]=-(-cos(angy)*a2 +sin(angx)*b1)/cos(angx+angy);

    L3[0]=-( sin(angy)*a2 -cos(angx)*b2)/cos(angx+angy);
    L3[1]=-(-cos(angy)*a2 +sin(angx)*b2)/cos(angx+angy);

    L4[0]=-( sin(angy)*a1 -cos(angx)*b2)/cos(angx+angy);
    L4[1]=-(-cos(angy)*a1 +sin(angx)*b2)/cos(angx+angy);


    /*Building the structured mesh*/
    Array<Mesh::Block> blks(1);
    blks[0].Set (/*NDim*/2, /*Tag*/-1, /*NVert*/4,
                 -1.,  L1[0], L1[1],  // tag, x, y, z
                 -2.,  L4[0], L4[1], 
                 -3.,  L3[0], L3[1], 
                 -4.,  L2[0], L2[1],
                 -10.,-20.,-30.,-40.); // face tags
    blks[0].SetNx (Nx);
    blks[0].SetNy (Ny);

    Mesh::Structured mesh(/*NDim*/2);
    mesh.Generate (blks,false);

    mesh.WriteVTU ("mymesh01");
    cout << " File <mymesh01.vtu> generated\n";

    /*DEM::Domain d01;
    d01.GenFromMesh(mesh,0.01,3.0,true,false,1.0);
    d01.WriteXDMF("dem_fracking_01");*/
	
    /*Checking for elements to move/delete*/
    /*Circling over the cells*/
    int i = 0;
    /*Taging cells*/
    while (i < mesh.Cells.Size())
    {
		mesh.Cells[i]->Tag = 0;
		for (int j=0; j<4; j++)
		{
			double m1 = mesh.Cells[i]->V[j]->C[0];
			double m2 = mesh.Cells[i]->V[j]->C[1];
			if (m1<0.0+lx/1000.0*cos(angx))
			{
				mesh.Cells[i]->V[j]->Tag=-40;
				mesh.Cells[i]->Tag++;
			}
			else if (m1>Lx-lx/1000.0*cos(angx))
			{
				mesh.Cells[i]->V[j]->Tag=-20;
				mesh.Cells[i]->Tag++;
			}
			else if (m2<0.0+ly/1000.0*cos(angy))
			{
				mesh.Cells[i]->V[j]->Tag=-10;
				mesh.Cells[i]->Tag++;
			}
			else if (m2>Ly-ly/1000.0*cos(angy))
			{
				mesh.Cells[i]->V[j]->Tag=-30;
				mesh.Cells[i]->Tag++;
			}
		}
    	i++;
    }
    
    i=0;
    while (i < mesh.Verts.Size())
    {
        /*Removing outside vertices*/
    	int fs=9;
    	mesh.Verts[i]->ID=i;
    	for (size_t j=0; j<mesh.Verts[i]->Shares.Size(); j++)
    	{
    		fs=min(fs,mesh.Verts[i]->Shares[j].C->Tag);
		}
    	if (fs > 3)
    	{
    		mesh.Verts.DelItem(i);
    		i--;
    	}
    	i++;
    }

    i=0;
    while (i < mesh.Cells.Size())
    {
    	mesh.Cells[i]->ID=i;
        /*Removing outside cells*/
    	if (mesh.Cells[i]->Tag == 4)
    	{
			Mesh::Cell * Cc = mesh.Cells[i];    		
			for (int j=0; j < 4; j++)
    		{
				int k=0;
				while (k < mesh.Cells[i]->V[j]->Shares.Size())
				{
					if (Cc==mesh.Cells[i]->V[j]->Shares[k].C)
					{
						mesh.Cells[i]->V[j]->Shares[k].C=0;
						mesh.Cells[i]->V[j]->Shares.DelItem(k);
						k--;
					}
					k++;
				}
			}
    		mesh.Cells.DelItem(i);
    		i--;
    	}

    	i++;
    }
    
    mesh.WriteVTU ("mymesh02");
    cout << " File <mymesh02.vtu> generated\n";
    
    /*DEM::Domain d02;
    d02.GenFromMesh(mesh,0.01,3.0,true,false,1.0);
    d02.WriteXDMF("dem_fracking_02");*/
    
    /*Fixing corner elements*/
    i=0;
    int ar1[4] = {-10, -10, -10, -10 } ;
    double ar2[4] = {sqrt(pow(Lx,2)+pow(Ly,2)), sqrt(pow(Lx,2)+pow(Ly,2)), 
					 sqrt(pow(Lx,2)+pow(Ly,2)), sqrt(pow(Lx,2)+pow(Ly,2))};
    double ar3[4] = {0.0+lx/100.0, Lx-lx/100.0, Lx-lx/100.0, 0.0+lx/100.0};
    double ar4[4] = {0.0+ly/100.0, 0.0-ly/100.0, Ly-ly/100.0, Ly-ly/100.0};
    double cc[4];
    double xx,yy;
	
    i=0;
    while (i < mesh.Cells.Size())
    {
		xx = (mesh.Cells[i]->V[0]->C[0]+mesh.Cells[i]->V[2]->C[0])/2.0;
		yy = (mesh.Cells[i]->V[0]->C[1]+mesh.Cells[i]->V[2]->C[1])/2.0;
		for (int j=0; j<4; j++)
		{
			cc[j] = sqrt(pow(mesh.Cells[i]->V[j]->C[0]-ar3[j],2)+pow(mesh.Cells[i]->V[j]->C[1]-ar4[j],2));
			if ((cc[j]<ar2[j])&&(mesh.Cells[i]->Tag<3))
			{
				ar2[j]=cc[j];
				ar1[j]=i;
			}
		}
		i++;
	}
	
	mesh.Cells[ar1[0]]->Tag  = -10;
	mesh.Cells[ar1[1]]->Tag  = -10;
	mesh.Cells[ar1[2]]->Tag  = -10;
	mesh.Cells[ar1[3]]->Tag  = -10;
	
	/*Handling corner (0,0)*/
	size_t k = 0;
	int m1 = 1;
	int m2 = 3;
	int m3 = 2;
	
	mesh.Cells[ar1[0]]->V[k]->C[0]=0.0;
	mesh.Cells[ar1[0]]->V[k]->C[1]=0.0;
	
    mesh.Cells[ar1[0]]->V[m3]->C[0]=max(mesh.Cells[ar1[0]]->V[m3]->C[0],lx/4.0);
    mesh.Cells[ar1[0]]->V[m3]->C[1]=max(mesh.Cells[ar1[0]]->V[m3]->C[1],ly/4.0);	
    
    a1=cos(angy)*mesh.Cells[ar1[0]]->V[m3]->C[0]+sin(angy)*mesh.Cells[ar1[0]]->V[m3]->C[1];
    mesh.Cells[ar1[0]]->V[m1]->C[0]=min(max(a1/cos(angy),lx/4.0),lx*cos(angx));
    mesh.Cells[ar1[0]]->V[m1]->C[1]=0.0;	    	    		
    
    mesh.Cells[ar1[0]]->V[m2]->C[0]=0.0;
    a1=sin(angx)*mesh.Cells[ar1[0]]->V[m3]->C[0]+cos(angx)*mesh.Cells[ar1[0]]->V[m3]->C[1];
    mesh.Cells[ar1[0]]->V[m2]->C[1]=min(max(a1/cos(angx),ly/4.0),ly*cos(angy));

    
	/*Marking edge nodes*/
	mesh.Cells[ar1[0]]->V[k]->Tag  = 0;
    
    int kk=mesh.Cells[ar1[0]]->V[k]->ID;
    for (size_t j=0; j<mesh.Verts[kk]->Shares.Size(); j++)
    {
		if (mesh.Verts[kk]->Shares[j].C->ID!=ar1[0])
		{
			mesh.Cells[mesh.Verts[kk]->Shares[j].C->ID]->Tag=4;
		}
	}
	
    kk=mesh.Cells[ar1[0]]->V[m1]->ID;
    for (size_t j=0; j<mesh.Verts[kk]->Shares.Size(); j++)
    {
		if ((mesh.Verts[kk]->Shares[j].C->Tag!=0)&&(mesh.Verts[kk]->Tag==0))
		{
			mesh.Verts[kk]->Shares[j].C->Tag++;
		}
	}	
	
	kk=mesh.Cells[ar1[0]]->V[m2]->ID;
    for (size_t j=0; j<mesh.Verts[kk]->Shares.Size(); j++)
    {
		if ((mesh.Verts[kk]->Shares[j].C->Tag!=0)&&(mesh.Verts[kk]->Tag==0))
		{
			mesh.Verts[kk]->Shares[j].C->Tag++;
		}
	}	
	
	mesh.Cells[ar1[0]]->V[m1]->Tag = -10;
    mesh.Cells[ar1[0]]->V[m2]->Tag = -40;
	
    /*Handling corner (Lx,0)*/
	k = 1;
	m1 = 2;
	m2 = 0;
	m3 = 3;
	
	mesh.Cells[ar1[1]]->V[k]->C[0]=Lx;
	mesh.Cells[ar1[1]]->V[k]->C[1]=0.0;
	
	mesh.Cells[ar1[1]]->V[m3]->C[0]=min(mesh.Cells[ar1[1]]->V[m3]->C[0],Lx-lx/4.0);
    mesh.Cells[ar1[1]]->V[m3]->C[1]=max(mesh.Cells[ar1[1]]->V[m3]->C[1],ly/4.0);	
    
    mesh.Cells[ar1[1]]->V[m1]->C[0]=Lx;
   	a1=sin(angx)*mesh.Cells[ar1[1]]->V[m3]->C[0]+cos(angx)*mesh.Cells[ar1[1]]->V[m3]->C[1];
    mesh.Cells[ar1[1]]->V[m1]->C[1]=min(max((a1-sin(angx)*Lx)/cos(angx),ly/4.0),ly*cos(angy));
    	    				
    a1=cos(angx)*mesh.Cells[ar1[1]]->V[m3]->C[0]+sin(angx)*mesh.Cells[ar1[1]]->V[m3]->C[1];
	mesh.Cells[ar1[1]]->V[m2]->C[0]=max(min(a1/cos(angx),Lx-lx/4.0),Lx-lx*cos(angx));
	mesh.Cells[ar1[1]]->V[m2]->C[1]=0.0;
	
	/*Marking edge nodes*/
	mesh.Cells[ar1[1]]->V[k]->Tag  = 0;

	
	
	kk=mesh.Cells[ar1[1]]->V[k]->ID;
    for (size_t j=0; j<mesh.Verts[kk]->Shares.Size(); j++)
    {
		if (mesh.Verts[kk]->Shares[j].C->ID!=ar1[1])
		{
			mesh.Cells[mesh.Verts[kk]->Shares[j].C->ID]->Tag=4;
		}
	}	
	
	kk=mesh.Cells[ar1[1]]->V[m1]->ID;
    for (size_t j=0; j<mesh.Verts[kk]->Shares.Size(); j++)
    {
		if ((mesh.Verts[kk]->Shares[j].C->Tag!=0)&&(mesh.Verts[kk]->Tag==0))
		{
			mesh.Verts[kk]->Shares[j].C->Tag++;
		}
	}	
	
	kk=mesh.Cells[ar1[1]]->V[m2]->ID;
    for (size_t j=0; j<mesh.Verts[kk]->Shares.Size(); j++)
    {
		if ((mesh.Verts[kk]->Shares[j].C->Tag!=0)&&(mesh.Verts[kk]->Tag==0))
		{
			mesh.Verts[kk]->Shares[j].C->Tag++;
		}
	}
	
    mesh.Cells[ar1[1]]->V[m1]->Tag = -20;
    mesh.Cells[ar1[1]]->V[m2]->Tag = -10;	
			
    /*Handling corner (Lx,Ly)*/
	k = 2;
	m1 = 3;
	m2 = 1;
	m3 = 0;
	
	mesh.Cells[ar1[2]]->V[k]->C[0]=Lx;
	mesh.Cells[ar1[2]]->V[k]->C[1]=Ly;
	
    mesh.Cells[ar1[2]]->V[m3]->C[0]=min(mesh.Cells[ar1[2]]->V[m3]->C[0],Lx-lx/4.0);
    mesh.Cells[ar1[2]]->V[m3]->C[1]=min(mesh.Cells[ar1[2]]->V[m3]->C[1],Ly-ly/4.0);	
    
    a1=cos(angy)*mesh.Cells[ar1[2]]->V[m3]->C[0]+sin(angy)*mesh.Cells[ar1[2]]->V[m3]->C[1];
   	mesh.Cells[ar1[2]]->V[m1]->C[0]=max(min((a1-sin(angy)*Ly)/cos(angy),Lx-lx/4.0),Lx-lx*cos(angx));
    mesh.Cells[ar1[2]]->V[m1]->C[1]=Ly; 
    	    				   	    				
    mesh.Cells[ar1[2]]->V[m2]->C[0]=Lx;
    a1=sin(angx)*mesh.Cells[ar1[2]]->V[m3]->C[0]+cos(angx)*mesh.Cells[ar1[2]]->V[m3]->C[1];
    mesh.Cells[ar1[2]]->V[m2]->C[1]=max(min((a1-sin(angx)*Lx)/cos(angx),Ly-ly/4.0),Ly-ly*cos(angy));
    
    /*Marking edge nodes*/
	mesh.Cells[ar1[2]]->V[k]->Tag  = 0;

    
    kk=mesh.Cells[ar1[2]]->V[k]->ID;
    for (size_t j=0; j<mesh.Verts[kk]->Shares.Size(); j++)
    {
		if (mesh.Verts[kk]->Shares[j].C->ID!=ar1[2])
		{
			mesh.Cells[mesh.Verts[kk]->Shares[j].C->ID]->Tag=4;
		}
	}
	
	kk=mesh.Cells[ar1[2]]->V[m1]->ID;
    for (size_t j=0; j<mesh.Verts[kk]->Shares.Size(); j++)
    {
		if ((mesh.Verts[kk]->Shares[j].C->Tag!=0)&&(mesh.Verts[kk]->Tag==0))
		{
			mesh.Verts[kk]->Shares[j].C->Tag++;
		}
	}	
	
	kk=mesh.Cells[ar1[2]]->V[m2]->ID;
    for (size_t j=0; j<mesh.Verts[kk]->Shares.Size(); j++)
    {
		if ((mesh.Verts[kk]->Shares[j].C->Tag!=0)&&(mesh.Verts[kk]->Tag==0))
		{
			mesh.Verts[kk]->Shares[j].C->Tag++;
		}
	}
	
    mesh.Cells[ar1[2]]->V[m1]->Tag = -30;
    mesh.Cells[ar1[2]]->V[m2]->Tag = -20;
    	    
    /*Handling corner (0,Ly)*/
	k = 3;
	m1 = 0;
	m2 = 2;
	m3 = 1;	
	
	mesh.Cells[ar1[3]]->V[k]->C[0]=0.0;
	mesh.Cells[ar1[3]]->V[k]->C[1]=Ly;
	
    mesh.Cells[ar1[3]]->V[m3]->C[0]=max(mesh.Cells[ar1[3]]->V[m3]->C[0],lx/4.0);
    mesh.Cells[ar1[3]]->V[m3]->C[1]=min(mesh.Cells[ar1[3]]->V[m3]->C[1],Ly-ly/4.0);	
    
	mesh.Cells[ar1[3]]->V[m1]->C[0]=0.0;
    a1=sin(angx)*mesh.Cells[ar1[3]]->V[m3]->C[0]+cos(angx)*mesh.Cells[ar1[3]]->V[m3]->C[1];
    mesh.Cells[ar1[3]]->V[m1]->C[1]=max(min(a1/cos(angx),Ly-ly/4.0),Ly-ly*cos(angy));
    	    				
    a1=cos(angx)*mesh.Cells[ar1[3]]->V[m3]->C[0]+sin(angx)*mesh.Cells[ar1[3]]->V[m3]->C[1];
    mesh.Cells[ar1[3]]->V[m2]->C[0]=min(max((a1-sin(angx)*Ly)/cos(angx),lx/4.0),lx*cos(angx));
    mesh.Cells[ar1[3]]->V[m2]->C[1]=Ly;
    
    /*Marking edge nodes*/
	mesh.Cells[ar1[3]]->V[k]->Tag  = 0;

	kk=mesh.Cells[ar1[3]]->V[k]->ID;
    for (size_t j=0; j<mesh.Verts[kk]->Shares.Size(); j++)
    {
		if (mesh.Verts[kk]->Shares[j].C->ID!=ar1[3])
		{
			mesh.Cells[mesh.Verts[kk]->Shares[j].C->ID]->Tag=4;
		}
	}    
    
     kk=mesh.Cells[ar1[3]]->V[m1]->ID;
    for (size_t j=0; j<mesh.Verts[kk]->Shares.Size(); j++)
    {
		if ((mesh.Verts[kk]->Shares[j].C->Tag!=0)&&(mesh.Verts[kk]->Tag==0))
		{
			mesh.Verts[kk]->Shares[j].C->Tag++;
		}
	}	
	
	kk=mesh.Cells[ar1[3]]->V[m2]->ID;
    for (size_t j=0; j<mesh.Verts[kk]->Shares.Size(); j++)
    {
		if ((mesh.Verts[kk]->Shares[j].C->Tag!=0)&&(mesh.Verts[kk]->Tag==0))
		{
			mesh.Verts[kk]->Shares[j].C->Tag++;
		}
	}
	
	mesh.Cells[ar1[3]]->V[m1]->Tag = -40;
    mesh.Cells[ar1[3]]->V[m2]->Tag = -30;
    
    
    /*Checking for elements to move/delete*/
    /*Circling over the cells*/
    i=0;
    while (i < mesh.Verts.Size())
    {
        /*Removing outside vertices*/
    	int fs=9;
    	mesh.Verts[i]->ID=i;
    	for (size_t j=0; j<mesh.Verts[i]->Shares.Size(); j++)
    	{
    		fs=min(fs,mesh.Verts[i]->Shares[j].C->Tag);
    	}
    	if (fs > 3)
    	{
    		mesh.Verts.DelItem(i);
    		i--;
    	}
    	i++;
    }

    i=0;
    while (i < mesh.Cells.Size())
    {
    	mesh.Cells[i]->ID=i;
        /*Removing outside cells*/
    	if (mesh.Cells[i]->Tag > 3)
    	{
			Mesh::Cell * Cc = mesh.Cells[i];    		
			for (int j=0; j < 4; j++)
    		{
				if (mesh.Cells[i]->V[j]->Tag==0)
				{
					double x40 = sqrt(pow(mesh.Cells[i]->V[j]->C[0],2));
					double x20 = sqrt(pow(mesh.Cells[i]->V[j]->C[0]-Lx,2));
					double x10 = sqrt(pow(mesh.Cells[i]->V[j]->C[1],2));
					double x30 = sqrt(pow(mesh.Cells[i]->V[j]->C[1]-Ly,2));
					double xx=x40;
					mesh.Cells[i]->V[j]->Tag=-40;
					if (xx>x20)
					{
						mesh.Cells[i]->V[j]->Tag=-20;
						xx=x20;
					}
					if (xx>x10)
					{
						mesh.Cells[i]->V[j]->Tag=-10;
						xx=x10;
					}
					if (xx>x30)
					{
						mesh.Cells[i]->V[j]->Tag=-30;
						xx=x30;
					}
					int k=0;
					while (k < mesh.Cells[i]->V[j]->Shares.Size())
					{
						mesh.Cells[i]->V[j]->Shares[k].C->Tag++;
						k++;
					}
				}
				int k=0;
				while (k < mesh.Cells[i]->V[j]->Shares.Size())
				{
					if (Cc==mesh.Cells[i]->V[j]->Shares[k].C)
					{
						mesh.Cells[i]->V[j]->Shares[k].C=0;
						mesh.Cells[i]->V[j]->Shares.DelItem(k);
						k--;
					}
					k++;
				}
			}
    		mesh.Cells.DelItem(i);

    		i--;
    	}

    	i++;
    }
    
    
    
    mesh.WriteVTU ("mymesh03");
    cout << " File <mymesh03.vtu> generated\n";
    
    /*DEM::Domain d03;
    d03.GenFromMesh(mesh,0.001,3.0,true,false,1.0);
    d03.WriteXDMF("dem_fracking_03"); */
    
  
    i=0;
    while (i < mesh.Cells.Size())
    {
        /*Correcting cells with tag 2*/
    	if (mesh.Cells[i]->Tag == 2)
    	{
    		/*Getting the coordinates of the vertices*/
    		for (size_t k=0; k<4; k++)
    		{
    			int ij=mesh.Cells[i]->V[k]->ID;
    			if (mesh.Verts[ij]->Tag!=0)
    			{
					int m1=k-1;
					int m2=k+1;    		    
					if (m1<0) 
						m1=3;
					else if (m2>3)
						m2=0;
					
					double cd[8];
    	            double cro[2];
    	            if (mesh.Cells[i]->V[m1]->Tag ==0)
    	            { 	    
						cd[0] = mesh.Cells[i]->V[m1]->C[0];    	    
						cd[4] = mesh.Cells[i]->V[m1]->C[1];
					}
					else
					{
						cd[0] = mesh.Cells[i]->V[m2]->C[0];    	    
						cd[4] = mesh.Cells[i]->V[m2]->C[1];
					}
    	    
					cd[1] = mesh.Cells[i]->V[k]->C[0];
					cd[5] = mesh.Cells[i]->V[k]->C[1];
   		    			
					if (mesh.Verts[ij]->Tag==-40)
					{
						cd[2] = 0.0;
						cd[3] = 0.0;
						cd[6] = 0.0;
						cd[7] = Ly;
					}

					if (mesh.Verts[ij]->Tag==-20)
					{
						cd[2] = Lx;
						cd[3] = Lx;
						cd[6] = 0.0;
						cd[7] = Ly;
					}
					if (mesh.Verts[ij]->Tag==-10)
					{
						cd[2] = 0.0;
						cd[3] = Lx;
						cd[6] = 0.0;
						cd[7] = 0.0;
					}
					if (mesh.Verts[ij]->Tag==-30)
					{
						cd[2] = 0.0;
						cd[3] = Lx;
						cd[6] = Ly;
						cd[7] = Ly;
					}
					
					cross(cd,cro);
					mesh.Verts[ij]->C[0]=cro[0];
					mesh.Verts[ij]->C[1]=cro[1];
				}
			}
		}
    	i++;
    }
    
    
    /*Checking for elements to move/delete*/
    /*Circling over the cells*/
    i=0;
    while (i < mesh.Verts.Size())
    {
        /*Removing outside vertices*/
    	int fs=9;
    	mesh.Verts[i]->ID=i;
    	for (size_t j=0; j<mesh.Verts[i]->Shares.Size(); j++)
    	{
    		fs=min(fs,mesh.Verts[i]->Shares[j].C->Tag);
    	}
    	if (fs > 3)
    	{
    		mesh.Verts.DelItem(i);
    		i--;
    	}
    	i++;
    }

    i=0;
    while (i < mesh.Cells.Size())
    {
    	mesh.Cells[i]->ID=i;
        /*Removing outside cells*/
    	if (mesh.Cells[i]->Tag > 3)
    	{
			Mesh::Cell * Cc = mesh.Cells[i];    		
			for (int j=0; j < 4; j++)
    		{
				int k=0;
				while (k < mesh.Cells[i]->V[j]->Shares.Size())
				{
					if (Cc==mesh.Cells[i]->V[j]->Shares[k].C)
					{
						mesh.Cells[i]->V[j]->Shares[k].C=0;
						mesh.Cells[i]->V[j]->Shares.DelItem(k);
						k--;
					}
					k++;
				}
			}
    		mesh.Cells.DelItem(i);

    		i--;
    	}

    	i++;
    }    
    
    
    mesh.WriteVTU ("mymesh04");
    cout << " File <mymesh04.vtu> generated\n";
    
    /*DEM::Domain d04;
    d04.GenFromMesh(mesh,0.001,3.0,true,false,1.0);
    d04.WriteXDMF("dem_fracking_04"); */    
    

    i=0;
   
     while (i<mesh.Cells.Size())
    {
		if (mesh.Cells[i]->Tag==1)
		{
			double xx=10000.0;
			int fs=0;
			for (size_t k=0; k<4; k++)
			{
				if (mesh.Cells[i]->V[k]->Tag < 0)
				{
					fs=k;
				}
			}
			switch(mesh.Cells[i]->V[fs]->Tag)
			{
				case -10:
				    xx = sqrt(pow(mesh.Cells[i]->V[fs]->C[1],2));
					break;
				case -20:
					xx = sqrt(pow(mesh.Cells[i]->V[fs]->C[0]-Lx,2));
					break;
				case -30:
					xx = sqrt(pow(mesh.Cells[i]->V[fs]->C[1]-Ly,2));
					break;
				case -40:
					xx = sqrt(pow(mesh.Cells[i]->V[fs]->C[0],2));
					break;			
			}
			if (xx<0.00001)
			{
				mesh.Cells[i]->Tag=0;
			}
		}
		i++;
	}
	
	i=0;
    while (i<mesh.Cells.Size())
    {
		if (mesh.Cells[i]->Tag==1)
		{
			int fs=0;
			/*Locating the outside vertex*/
			for (size_t k=0; k<4; k++)
			{
				if (mesh.Cells[i]->V[k]->Tag < 0)
				{
					fs=k;
				}
			}
			int ij=mesh.Cells[i]->V[fs]->ID;
			mesh.Cells[i]->V[fs]->Tag=10;
			
			/*Locating the inside vertex neigbors*/
			int m1=fs-1;
			int m2=fs+1;
			if (m1<0)
				m1=3;
			else if (m2>3)
				m2=0;
			int m3 = (fs+2+4)%4;	
				
				
			/* Duplicating outside Vertex*/
			Mesh::Vertex * v      = new Mesh::Vertex;
			v->ID=mesh.Verts.Size();
			v->Tag = 5;
			v->C[0] = mesh.Verts[ij]->C[0];
			v->C[1] = mesh.Verts[ij]->C[1];
			int iv = v->ID;
			mesh.Verts.Push(v);
			
			/* Duplicating the Cell*/
			Mesh::Cell * e      = new Mesh::Cell;
			e->ID=mesh.Cells.Size();
			e->Tag=0;
			int ie = e->ID;
			e->V.Resize(4);
			mesh.Cells.Push(e);
			
			/* Duplicating the Cell*/
			Mesh::Cell * a      = new Mesh::Cell;
			a->ID=mesh.Cells.Size();
			a->Tag=0;
			int ia = a->ID;
			a->V.Resize(4);
			mesh.Cells.Push(a);
			
			double cd[8];
    	    double cro[2];
    	    
    		/*Changing the Vertex coordinates*/
			if (mesh.Verts[ij]->C[0]<0.0+0.0)
			{
				cd[0] = mesh.Cells[i]->V[fs]->C[0];
				cd[4] = mesh.Cells[i]->V[fs]->C[1];
			    cd[2] = 0.0;
			    cd[3] = 0.0;
    	        cd[6] = 0.0;
    	        cd[7] = Ly;
    	        
    	        cd[1] = mesh.Cells[i]->V[m1]->C[0];
                cd[5] = mesh.Cells[i]->V[m1]->C[1];
                cross(cd,cro);
                mesh.Verts[ij]->C[0]=cro[0];
                mesh.Verts[ij]->C[1]=cro[1];
                
                cd[1] = mesh.Cells[i]->V[m2]->C[0];
                cd[5] = mesh.Cells[i]->V[m2]->C[1];
                cross(cd,cro);
                mesh.Verts[iv]->C[0]=cro[0];
                mesh.Verts[iv]->C[1]=cro[1];
                
                if (mesh.Cells[i]->V[m1]->C[0]<0.0+lx/4*cos(angx))
                {
					mesh.Cells[i]->V[m1]->C[0]=0.0+lx/4*cos(angx);
				}
				if (mesh.Cells[i]->V[m2]->C[0]<0.0+lx/4*cos(angx))
                {
					mesh.Cells[i]->V[m2]->C[0]=0.0+lx/4*cos(angx);
				}
			}
			
			if (mesh.Verts[ij]->C[0]>Lx-0.0)
			{
           	    cd[0] = mesh.Cells[i]->V[fs]->C[0];
				cd[4] = mesh.Cells[i]->V[fs]->C[1];
                cd[2] = Lx;
			    cd[3] = Lx;
    	        cd[6] = 0.0;
    	        cd[7] = Ly;
    	        
    	        cd[1] = mesh.Cells[i]->V[m1]->C[0];
                cd[5] = mesh.Cells[i]->V[m1]->C[1];
                cross(cd,cro);
                mesh.Verts[ij]->C[0]=cro[0];
                mesh.Verts[ij]->C[1]=cro[1];
                
                cd[1] = mesh.Cells[i]->V[m2]->C[0];
                cd[5] = mesh.Cells[i]->V[m2]->C[1];
                cross(cd,cro);
                mesh.Verts[iv]->C[0]=cro[0];
                mesh.Verts[iv]->C[1]=cro[1];
                
                if (mesh.Cells[i]->V[m1]->C[0]>Lx-lx/4*cos(angx))
                {
					mesh.Cells[i]->V[m1]->C[0]=Lx-lx/4*cos(angx);
				}
				if (mesh.Cells[i]->V[m2]->C[0]>Lx-lx/4*cos(angx))
                {
					mesh.Cells[i]->V[m2]->C[0]=Lx-lx/4*cos(angx);
				}
			}
			if (mesh.Verts[ij]->C[1]<0.0+0.0)
			{
				cd[0] = mesh.Cells[i]->V[fs]->C[0];
    	        cd[4] = mesh.Cells[i]->V[fs]->C[1];
				cd[2] = 0.0;
			    cd[3] = Lx;
    	        cd[6] = 0.0;
    	        cd[7] = 0.0;
    	        
    	        cd[1] = mesh.Cells[i]->V[m1]->C[0];
                cd[5] = mesh.Cells[i]->V[m1]->C[1];
                cross(cd,cro);
                mesh.Verts[ij]->C[0]=cro[0];
                mesh.Verts[ij]->C[1]=cro[1];
                
                cd[1] = mesh.Cells[i]->V[m2]->C[0];
                cd[5] = mesh.Cells[i]->V[m2]->C[1];
                cross(cd,cro);
                mesh.Verts[iv]->C[0]=cro[0];
                mesh.Verts[iv]->C[1]=cro[1];
                
                if (mesh.Cells[i]->V[m1]->C[1]<0.0+ly/4*cos(angy))
                {
					mesh.Cells[i]->V[m1]->C[1]=0.0+ly/4*cos(angy);
				}
				if (mesh.Cells[i]->V[m2]->C[1]<0.0+ly/4*cos(angy))
                {
					mesh.Cells[i]->V[m2]->C[1]=0.0+ly/4*cos(angy);
				}
			}
			if (mesh.Verts[ij]->C[1]>Ly-0.0)
			{
                cd[0] = mesh.Cells[i]->V[fs]->C[0];
				cd[4] = mesh.Cells[i]->V[fs]->C[1];    
                cd[2] = 0.0;
			    cd[3] = Lx;
    	        cd[6] = Ly;
    	        cd[7] = Ly;
    	        
    	        cd[1] = mesh.Cells[i]->V[m1]->C[0];
                cd[5] = mesh.Cells[i]->V[m1]->C[1];
                cross(cd,cro);
                mesh.Verts[ij]->C[0]=cro[0];
                mesh.Verts[ij]->C[1]=cro[1];
                
                cd[1] = mesh.Cells[i]->V[m2]->C[0];
                cd[5] = mesh.Cells[i]->V[m2]->C[1];
                cross(cd,cro);
                mesh.Verts[iv]->C[0]=cro[0];
                mesh.Verts[iv]->C[1]=cro[1];
                
                if (mesh.Cells[i]->V[m1]->C[1]>Ly-ly/4*cos(angy))
                {
					mesh.Cells[i]->V[m1]->C[1]=Ly-ly/4*cos(angy);
				}
				if (mesh.Cells[i]->V[m2]->C[1]>Ly-ly/4*cos(angy))
                {
					mesh.Cells[i]->V[m2]->C[1]=Ly-ly/4*cos(angy);
				}
			}
					
			
			/* Generating middle Vertex*/
			Mesh::Vertex * vm      = new Mesh::Vertex;
			vm->ID=mesh.Verts.Size();
			vm->Tag=2;
			int im = vm->ID;
			mesh.Verts.Push(vm);			

			mesh.Verts[im]->C[0]=(mesh.Verts[ij]->C[0]+mesh.Verts[iv]->C[0])/2.0;
			mesh.Verts[im]->C[1]=(mesh.Verts[ij]->C[1]+mesh.Verts[iv]->C[1])/2.0;

            
			/*Sorting vertices to cells*/
			/*ij goes with m1 and belongs to i*/
			/*Substituting im for m2 in elements containing m1*/
			/*int ls = mesh.Cells[i]->V[m2]->ID;
			mesh.Cells[i]->V[m2]=0;
			mesh.Cells[i]->V[m2]=mesh.Verts[im];
			mesh.Cells[i]->V[m2]=mesh.Verts[iv];
			mesh.Cells[i]->V[m2]=mesh.Verts[ls];*/
			

			mesh.Cells[ie]->V[m3]=mesh.Cells[i]->V[m3];
			mesh.Cells[ie]->V[fs]=mesh.Verts[ij];
			mesh.Cells[ie]->V[m1]=mesh.Cells[i]->V[m1];
			mesh.Cells[ie]->V[m2]=mesh.Verts[im];
			
			int k=0;
			while (k < mesh.Cells[i]->V[m1]->Shares.Size())
			{
				if (mesh.Cells[i]==mesh.Cells[i]->V[m1]->Shares[k].C)
				{
					mesh.Cells[i]->V[m1]->Shares[k].C=mesh.Cells[ie];
				}
				k++;
			}


			mesh.Cells[ia]->V[m3]=mesh.Cells[i]->V[m3];
			mesh.Cells[ia]->V[fs]=mesh.Verts[iv];
			mesh.Cells[ia]->V[m1]=mesh.Verts[im];
			mesh.Cells[ia]->V[m2]=mesh.Cells[i]->V[m2];
			
			k=0;
			while (k < mesh.Cells[i]->V[m2]->Shares.Size())
			{
				if (mesh.Cells[i]==mesh.Cells[i]->V[m2]->Shares[k].C)
				{
					mesh.Cells[i]->V[m2]->Shares[k].C=mesh.Cells[ia];
				}
				else
				{
					for (int j=0; j<4; j++)
					{
						if (mesh.Cells[i]->V[m2]->Shares[k].C->V[j]==mesh.Verts[ij])
							mesh.Cells[i]->V[m2]->Shares[k].C->V[j]=mesh.Verts[iv];
					}
				}
				k++;
			}

			mesh.Cells[i]->Tag=4;
			
		}
		i++;
	}
	
	

	i=0;
    while (i < mesh.Verts.Size())
    {
    	if (mesh.Verts[i]->Tag == 0)
    	{
			if (mesh.Verts[i]->C[0]<lx/8.0*cos(angx))
			{
				mesh.Verts[i]->C[0]=lx/8.0*cos(angx);
			}
			if (mesh.Verts[i]->C[0]>Lx-lx/8.0*cos(angx))
			{
				mesh.Verts[i]->C[0]=Lx-lx/8.0*cos(angx);
			}
			if (mesh.Verts[i]->C[1]<ly/8.0*cos(angy))
			{
				mesh.Verts[i]->C[1]=ly/8.0*cos(angy);
			}
			if (mesh.Verts[i]->C[1]>Ly-ly/8.0*cos(angy))
			{
				mesh.Verts[i]->C[1]=Ly-ly/8.0*cos(angy);
			}
    	}

    	i++;
    }
    
    i=0;
    while (i < mesh.Cells.Size())
    {
    	mesh.Cells[i]->ID=i;
        /*Removing outside cells*/
    	if (mesh.Cells[i]->Tag == 4)
    	{
			Mesh::Cell * Cc = mesh.Cells[i];    		
			for (int j=0; j < 4; j++)
    		{
				int k=0;
				while (k < mesh.Cells[i]->V[j]->Shares.Size())
				{
					if (Cc==mesh.Cells[i]->V[j]->Shares[k].C)
					{
						mesh.Cells[i]->V[j]->Shares[k].C=0;
						mesh.Cells[i]->V[j]->Shares.DelItem(k);
						k--;
					}
					k++;
				}
			}
    		mesh.Cells.DelItem(i);
    		i--;
    	}

    	i++;
    }

    i=0;
    /*checking for convexity*/
    while (i < mesh.Cells.Size())
    {
    	/*finding the intersection point of the diagonals*/
    	double cd[8];
    	double cro[2];
        cd[0] = mesh.Cells[i]->V[0]->C[0];
        cd[1] = mesh.Cells[i]->V[2]->C[0];
        cd[2] = mesh.Cells[i]->V[1]->C[0];
        cd[3] = mesh.Cells[i]->V[3]->C[0];
        cd[4] = mesh.Cells[i]->V[0]->C[1];
        cd[5] = mesh.Cells[i]->V[2]->C[1];
        cd[6] = mesh.Cells[i]->V[1]->C[1];
        cd[7] = mesh.Cells[i]->V[3]->C[1];
        cross(cd,cro);
    	double xx = cro[0];
    	double yy = cro[1];

   		double c1[9];
       	double ar1[4];

        c1[0]=xx;
        c1[1]=yy;
        c1[2]=1.0;
        c1[5]=1.0;
        c1[8]=1.0;

        for (size_t k=0; k<4; k++)
        {
        	c1[3]=mesh.Cells[i]->V[k]->C[0];
        	c1[4]=mesh.Cells[i]->V[k]->C[1];
        	c1[6]=mesh.Cells[i]->V[(k+1)%4]->C[0];
        	c1[7]=mesh.Cells[i]->V[(k+1)%4]->C[1];

        	ar1[k]=determ(c1);
        }
        if (ar1[0]<0 && ar1[1]<0 && ar1[2]<0 && ar1[3]<0)
        {
    		/*cout <<i<<" ! negative!"<<"\n";*/
        }
        else if (ar1[0]>0 && ar1[1]>0 && ar1[2]>0 && ar1[3]>0)
        {
    		/*cout <<i<<" ! positive!"<<"\n";*/
        }
        else
        {
    		cout <<i<<" ! non convex cell! Tag= "<<mesh.Cells[i]->Tag<<"\n";
    	}
    	i++;
    }
    cout <<"Cells with Tag=3 are to be modified now"<<"\n";

    mesh.WriteVTU ("mymesh05");
    cout << " File <mymesh05.vtu> generated\n";
    
    /*DEM::Domain d05;
    d05.GenFromMesh(mesh,0.001,3.0,true,false,1.0);
    d05.WriteXDMF("dem_fracking_05");*/
	
    cout << "number of Cells after 5: "<< mesh.Cells.Size() << "\n";
	i=0;
    while (i < mesh.Cells.Size())
    {
    	mesh.Cells[i]->ID=i;
        /*Removing outside cells*/
    	if (mesh.Cells[i]->Tag == 3)
    	{  		
			for (int j=0; j < 4; j++)
    		{
				if (mesh.Cells[i]->V[j]->Tag==0)
				{
					int m1=(j-1+4)%4;
					int m2=(j+1+4)%4;
					int m3=(j+2+4)%4;

					double cd[8];
					double cro[2];
					cd[0] = mesh.Cells[i]->V[m3]->C[0];
					cd[4] = mesh.Cells[i]->V[m3]->C[1];

					cd[1] = mesh.Cells[i]->V[j]->C[0];
					cd[5] = mesh.Cells[i]->V[j]->C[1];
    	    
					/*Changing the Vertex coordinates*/
					if (mesh.Cells[i]->V[m3]->Tag==-40)
					{
						cd[2] = 0.0;
						cd[3] = 0.0;
						cd[6] = 0.0;
						cd[7] = Ly;
					}
			
					if (mesh.Cells[i]->V[m3]->Tag==-20)
					{
						cd[2] = Lx;
						cd[3] = Lx;
						cd[6] = 0.0;
						cd[7] = Ly;
					}
					if (mesh.Cells[i]->V[m3]->Tag==-10)
					{
						cd[2] = 0.0;
						cd[3] = Lx;
						cd[6] = 0.0;
						cd[7] = 0.0;
					}
					if (mesh.Cells[i]->V[m3]->Tag==-30)
					{
						cd[2] = 0.0;
						cd[3] = Lx;
						cd[6] = Ly;
						cd[7] = Ly;
					}			
					                     
					cross(cd,cro);
					
					m1=mesh.Cells[i]->V[m1]->ID;
					m2=mesh.Cells[i]->V[m2]->ID;
					m3=mesh.Cells[i]->V[m3]->ID;
					
					mesh.Verts[m3]->Tag=-100;

					mesh.Verts[m1]->C[0]=cro[0];
					mesh.Verts[m1]->C[1]=cro[1];
					mesh.Verts[m2]->C[0]=cro[0];
					mesh.Verts[m2]->C[1]=cro[1];

					/*int k=0;
					while (k < mesh.Cells.Size())
					{ 
						for (int m=0; m < 4; m++)
						{					
							if ((mesh.Cells[k]->V[m]==mesh.Verts[m1])&&(mesh.Cells[k]->Tag!=3))
							{
								mesh.Cells[k]->V[m]=mesh.Cells[i]->V[m3];
								cout << k<<"one out"<<i<<" a "<<mesh.Verts[m1]->ID<<"\n";
							}
						}
						k++;
					}*/
				}
			}
    	}
    	i++;
    }
    
    /*i=0;
    while (i < mesh.Verts.Size())
    {
    	if (mesh.Verts[i]->Tag==-100)
    	{
    		mesh.Verts.DelItem(i);
    		i--;
    	}
    	i++;
    }*/
    
    i=0;
    while (i < mesh.Cells.Size())
    {
    	mesh.Cells[i]->ID=i;
        /*Removing outside cells*/
    	if (mesh.Cells[i]->Tag == 3)
    	{
			Mesh::Cell * Cc = mesh.Cells[i];    		
			for (int j=0; j < 4; j++)
    		{
				int k=0;
				while (k < mesh.Cells[i]->V[j]->Shares.Size())
				{
					if (Cc==mesh.Cells[i]->V[j]->Shares[k].C)
					{
						mesh.Cells[i]->V[j]->Shares[k].C=0;
						mesh.Cells[i]->V[j]->Shares.DelItem(k);
						k--;
					}
					k++;
				}
			}
    		mesh.Cells.DelItem(i);
    		i--;
    	}

    	i++;
    }
	
	
	i=0;
    while (i < mesh.Verts.Size())
    {
		if (mesh.Verts[i]->Tag == 0)
		{
			if (mesh.Verts[i]->Shares.Size()<3)
			{
				int k=0;
				for (int j=0; j < 4; j++) 
				{
					k=min(k,mesh.Verts[i]->Shares[0].C->V[j]->Tag);
				}
				
				mesh.Verts[i]->Tag = k;
				switch(mesh.Verts[i]->Tag)			
				{				
					case -10:
						mesh.Verts[i]->C[1]=0.0;
						break;
					case -20:
						mesh.Verts[i]->C[0]=Lx;
						break;				
					case -30:					
						mesh.Verts[i]->C[1]=Ly;
						break;
					case -40:
						mesh.Verts[i]->C[0]=0.0;
						break;			
				}
			}
		}
		i++;
	}
	
	
	
    mesh.WriteVTU ("mymesh06");
    cout << " File <mymesh06.vtu> generated\n";
    

    i=0;
    /*checking for convexity*/
    while (i < mesh.Cells.Size())
   	{
    	/*finding the intersection point of the diagonals*/
        double cd[8];
        double cro[2];

        cd[0] = mesh.Cells[i]->V[0]->C[0];
        cd[1] = mesh.Cells[i]->V[2]->C[0];
        cd[2] = mesh.Cells[i]->V[1]->C[0];
        cd[3] = mesh.Cells[i]->V[3]->C[0];
        cd[4] = mesh.Cells[i]->V[0]->C[1];
        cd[5] = mesh.Cells[i]->V[2]->C[1];
        cd[6] = mesh.Cells[i]->V[1]->C[1];
        cd[7] = mesh.Cells[i]->V[3]->C[1];
        cross(cd,cro);
        double xx = cro[0];
        double yy = cro[1];

        double c1[9];
        double ar1[4];

        c1[0]=xx;
        c1[1]=yy;
        c1[2]=1.0;
        c1[5]=1.0;
        c1[8]=1.0;

        for (size_t k=0; k<4; k++)
       	{
        	c1[3]=mesh.Cells[i]->V[k]->C[0];
          	c1[4]=mesh.Cells[i]->V[k]->C[1];
            c1[6]=mesh.Cells[i]->V[(k+1)%4]->C[0];
            c1[7]=mesh.Cells[i]->V[(k+1)%4]->C[1];

            ar1[k]=determ(c1);
         }
        if (ar1[0]<0 && ar1[1]<0 && ar1[2]<0 && ar1[3]<0)
        {
        	cout <<i<<" ! negative!"<<"\n";
        }
        else if (ar1[0]>0 && ar1[1]>0 && ar1[2]>0 && ar1[3]>0)
        {
        	/*cout <<i<<" ! positive!"<<"\n";*/
        }
        else
        {
        	cout <<i<<" ! non convex cell! Tag= "<<mesh.Cells[i]->Tag<<"\n";
        }
        i++;
    }		    
	
	/*Preparing to introduce dispersion*/
	i = 0;
    /*Taging the boundary cells*/
    while (i < mesh.Cells.Size())
    {
		mesh.Cells[i]->Tag = 0;
		for (int j=0; j<4; j++)
		{
			if (mesh.Cells[i]->V[j]->Tag!=0)
			{
				mesh.Cells[i]->Tag++;
			}
		}
    	i++;
    }
    
    i=0;
    /*Taging nodes belonging to boundary cells*/
    while (i < mesh.Verts.Size())
    {
		mesh.Verts[i]->Tag = 0;
		for (int j=0; j<mesh.Verts[i]->Shares.Size(); j++)
		{
			if (mesh.Verts[i]->Shares[j].C->Tag>0)
			{
				mesh.Verts[i]->Tag++;
			}
		}
    	i++;
    }
    
    /*Checking that the dispersion is small enough not to cause problems*/
    double deltax = min(sqrt(pow(cell1[0],2)), sqrt(pow(cell2[0],2)));
    deltax = min(deltax, sqrt(pow(cell3[0],2)));
    
    double deltay = min(sqrt(pow(cell1[1],2)), sqrt(pow(cell2[1],2)));
    deltay = min(deltay, sqrt(pow(cell3[1],2)));
    
    deltax =min(lx,ly);
    if (deltax < lsx*3.0)
    {
		lsx=0.0;
		cout<<"Value of x-dispersion is too big, setting it to 0 \n";
	}
	
	if (deltax < lsy*3.0)
    {
		lsy=0.0;
		cout<<"Value of y-dispersion is too big, setting it to 0 \n";
	}
	
	/*Applying the dispersion*/
	i=0;
	srand(time(0));
	double Big=static_cast<double>(RAND_MAX)/2.0;
    while (i < mesh.Verts.Size())
    {
		/*if (mesh.Verts[i]->Tag<1)*/
                 xx=sqrt(pow(mesh.Verts[i]->C[0]-Lx/2.0,2));
                 yy=sqrt(pow(mesh.Verts[i]->C[1]-Ly/2.0,2));
                 if (xx<2.0*Lx/6.0 && yy<2.0*Ly/6.0) 
		{
                  

			double r = (double)(rand()-Big)/Big;
			mesh.Verts[i]->C[0]=mesh.Verts[i]->C[0]+r*lsx;
		    r = (double)(rand()-Big)/Big;
			mesh.Verts[i]->C[1]=mesh.Verts[i]->C[1]+r*lsy;
		}
    	i++;
    }
    
    ofstream cells;
    cells.open ("cells.mesh");
    i=0;
    while (i < mesh.Cells.Size())
    {
		cells <<4<<setw(6)<<mesh.Cells[i]->V[0]->ID <<setw(6)<<mesh.Cells[i]->V[1]->ID <<
				   setw(6)<<mesh.Cells[i]->V[2]->ID <<setw(6)<<mesh.Cells[i]->V[3]->ID <<"\n";
		i++;
	}
    cells.close();
    
    ofstream verts;
    verts.open ("verts.mesh");
    
    i=0;
    while (i < mesh.Verts.Size())
    {
		verts <<mesh.Verts[i]->ID<<setw(12)<<mesh.Verts[i]->C[0] <<setw(12)<<mesh.Verts[i]->C[1] <<"\n";
		i++;
	}
    verts.close();
    
    mesh.WriteVTU ("mymesh00");
    mesh.WriteVTU ("mymesh07");
    cout << " File <mymesh07.vtu> generated\n";
}


MECHSYS_CATCH

